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ABSTRACT 

. We study the effect of disc self-gravity on instabilities associated with gaps opened 

by a giant Saturn mass planet in a protoplanetary disc that lead to the formation 
of vortices. We also study the nonlinear evolution of the vortices when this kind of 
instability occurs in a self-gravitating disc as well as the potential effect on type III 
planetary migration due to angular momentum exchange via coorbital flows. 

It is shown analytically and is confirmed through linear calculations that vortex 
£j . forming modes with low azimuthal mode number , m, are stabilised by the effect of self- 

gravity if the background structure is assumed fixed. However, the disc's self-gravity 
also affects the background gap surface density profile in a way that destabilises modes 
with high m. Linear calculations show that the combined effect of self-gravity through 
its effect on the background structure and its direct effect on the linear modes shifts 
the most rapidly growing vortex mode to higher m. 
lO ' Hydrodynamic simulations of gaps opened by a Saturn mass planet show more 

C^l \ vortices develop with increasing disc mass and therefore importance of self-gravity. 

For sufficiently large disc mass the vortex instability is suppressed, consistent with 
analytical expectations. In this case a new global instability develops instead. 
\ In the non-linear regime, we found that vortex merging is in general increasingly 

delayed as the disc mass increases and in some cases multiple vortices persist until 
the end of simulations. For massive discs in which the vortices merge, the post-merger 
vortex is localised in azimuth and has similar structure to a Kida-like vortex. This 
is unlike the case without self-gravity where vortices merge to form a larger vortex 
extended in azimuth. 

y\ . In order to study the properties of the vortex systems without the influence of the 

H ■ planet, we also performed a series of supplementary simulations of co-orbital Kida- 

like vortices. We found that self-gravity enables Kida-like vortices to execute horseshoe 
turns upon encountering each other. As a result vortex merging is avoided on time- 
scales where it would occur without self-gravity. Thus we suggest that mutual repulsion 
of self-gravitating vortices in a rotating flow is responsible for the delayed vortex 
merging seen in the disc-planet simulations. 

The effect of self-gravity on vortex-induced migration in low viscosity discs is 
briefly discussed. We found that when self-gravity is included and the disc mass is in the 
ra nge where vortex formin g instabilities occur, the vortex-induced type III migration 
of lLin fc Papaloizoul (|201Qh is delayed. There are also expected to be longer periods of 
slow migration between the short bursts of rapid migration compared to what occurs 
in a simulation without self-gravity. However, the extent of induced rapid migration is 
unchanged and involves flow of vortex material across the gap, independent of whether 
or not self-gravity is included. 
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1 INTRODUCTION 



Since the first detection of an extrasolar gian t planet orbit- 
i ng a m ain sequence star in a 4 day orbit bv lMavor &; Quelozl 
1995), there have been many observatio ns of hot Jupiters 



Marcv et ai1l2005l : lUdrv fc Santosll200/t ). One possible ex- 



planation for their close proximity to their host stars is or- 
bital migration due to interaction with the gaseous proto- 
planetary disc, fro m their site of formation to their current 
observed location. (Lin fc Papaloizoulll986l ) 

Planetary migr ati on was first studied by 
iGoldreich fe Tremaind ([l980), and since then many 
analytical and numerical stu dies of disc-planet inter actions 
have been undertaken (see Papalo izou et al.l 2007, for a 
review). However, most of these adopt a low mass disc 
model for which it may be assumed that self-gravity may 
be neglected. This is even the case when rapid type III 
migration that requires a massive disc is considered. It is 
the purpose of this paper to investigate the effect of the 
disc self-gravity on such disc planet interactions and in 
particular to consider its effect on the stability of the disc 
to vortex formation, a phenomenon that has been found to 
occur in discs with very low viscosity. 

Astrophysical discs may be dynamically unstable for 
a variety of reasons. One example is instability associated 
with steep surface density g r adien t s or vortenshrjj] extrema 
dPapaloizou fe Prinrid [l985l, Il987l: IPapaloizou fc^Linl [l989: 
[ Papaloizou &; Savoniid Il99ll ; lLovelace et al.lll999l : iLi et al.l 
2000 ) that can lead to vortex formation in non- linear regime 
(|Li et all l200ll ) . The disc vortensity profile necessary for 
such instabilities to occur may be induced by disc-planet 
interaction. For sufficiently large planetary masses and/or 
sufficiently low disc viscosity, the tidal perturbation on the 
gaseous disc leads to a dip or gap in the surface density 
([Lin fe Papaloizoul [l986). There have been a series of stud- 
i es of the vortex instability assoc i ated with planetary gaps 
([Koller et al.ll2003l : ILi et alJl2005l : Ide Val-Borro et al.ll2007l) 
and its consequences for planetary migration (|Ou et al.l 
120071 : ILi et al.ll2009l : iLin fe Papaloizoull2010l : lYu et alfeoTof T 
However, the effect of disc self- gravity is still unclear. 

In the context of disc-planet interactions, disc self- 
gravity is often neglected. It has thus far only been 
studie d i n a limited numb er of works (INelson &: Ben: 
2003allbl: IZhang et al.1 120081 : iBaruteau fe Massetl l20Qi 
Li et al. ( 2009) included self- gravity in their study of Type I 



migration in low viscosity discs, where the vortex instability 
develops for a viscosity parameters ^ 10 -5 , but the specific 
effect of self-gravity on the instability was not discussed an- 
alytically or explored numerically. In this paper we examine 
the effect of self-gravity on the vortex instability associated 
with edges of gaps opened by a Saturn-mass planet and its 
effect on the subsequent nonlinear evolution of vortices to- 
gether with a first calculation of the subsequent effects on 
orbital migration. 

This paper is organised as follows. In J2] we present 
the governing equations and models for our disc-planet sys- 
tems. We discuss the linearised stability problem in 33 We 
show that vortex forming modes with long azimuthal wave- 
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1 The term vortensity is used for the ratio of vorticity to surface 

density 



length associated with vortensity minima at gap edges, in- 
duced by an embedded planet, are stabilised by self- gravity 
when the background model is held fixed. In 3Hwe perform 
linear mode calculations for disc models of varying mass. 
We confirm the analytic discussion but show in addition 
that changes to the gap structure that occur with increas- 
ing disc mass are destabilising, an effect that is stronger 
for modes with short wavelength in the azimuthal direction. 
All of these effects cause the most unstable modes to shift 
to shorter azimuthal wavelength as the disc mass increases. 
However, in practice these modes are eventually dominated 
by long azimuthal wavelength, global edge modes, which are 
not associ ated with localised v ortex formation and which are 
studied in ILin & Papaloizou! ([201 ll ). Edge modes dominate 
our models when the disc-to-star mass ratio is > 0.047, cor- 
responding to a Keplerian Toomre stability parameter < 2 
at the outer disc boundary (see ^2.11 and Appendix [A]) . As 
this paper is concerned with vortex formation, we focus on 
models with disc mass below this threshold. 

We go on to present nonlinear hydrodynamic simula- 
tions in 33 In accord with expectations from linear theory 
the instability produces more vortices of smaller scale as the 
disc mass increases. The vortices are found to survive for 
much longer against merging than in the non self- gravitating 
limit and themselves become self- gravitating with strong 
over surface densities. This is found to enable pairs of vor- 
tices to undergo co-orbital dynamics with horseshoe trajec- 
tories which aids to inhibit merging. In 3S]we perform sim- 
ulations of interacting Kida-like vortices without an embed- 
ded planet in order to clarify results obtained for vortices 
at gap edges. These show that the co-orbital dynamics oc- 
curs independently of the embedded planet. In 33 we give a 
preliminary investigation of the effects of vortex formation 
in low viscosity discs on orbital migration when self-gravity 
is significant but not strong enough to prevent effective vor- 
tex formation. As in the non self- gravitating case, episodes of 
rapid migration occur as vortices are scattered by the planet 
([Lin fe Papaloizou! [2010h . but longer time intervals between 
them are expected on account of stabilisation by self-gravity 
and slowed vortex merging. Finally we summarise and con- 
clude in M 



2 DISC MODEL AND BASIC EQUATIONS 

We consider a gaseous disc of mass Ma orbiting a central star 
of mass M*. We assume the disc to be razor thin and thus 
work in the the two dimensional flat disc approximation. 
We adopt cylindrical polar co-ordinates (r, ip) centred on 
the star and defined in the plane of the disc. The reference 
frame is non-rotating. 

The governing hydrodynamic equations are the conti- 
nuity and momentum equations 



DT, 
~Dt 
Du 
~Dt 



-EV • u, 
1 



Vp - v$ + 



(1) 

(2) 



where £ is the surface density, u is the velocity field and p 
is the vertically integrated pressure. Numerical calculations 
adopt a locally isothermal equation of state, p = d?E where 
c s = HQk is the local sound speed and = \J GM* / r 3 = 
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2tt/P. Here H = hr is the putative disc semi-thickness. We 
fix the constant h = 0.05. The viscous force per unit area is / 
being characterised by a uniform kinematic viscosity v. The 
gravitational potential 4> includes the indirect potential 
the potential due to the central star, <I>* = —GM*/r, the 
potential due to the disc &d when self-gravity is included 
and the potential due to an embedded planet, < E > P , for disc- 
planet systems. Details of / are given in Masset (2002). The 
potentials due to the planet and the disc are 



GM V 



yjr 2 + r$- 2rr p cos (tp - ip p ) + e| 
GS(r',y') 



(3) 



L 



t> \/r 2 + r' 2 — 2rr' cos {Jp — <p*) + e 2 



r'dr'dcp', (4) 



where the integral is taken over the domain of the disc V. 
The planet has mass M p , position (r p , ip p ) and its potential 
is softened by use of the softening length e p = 0.6H(r p ). 
Similarly the disc potential is softened by use of the softening 
length e g = 0.3H(r'). The indirect potential takes account 
of the the forces on the central star, which is at the origin of 
a non inert ial frame, due to the disc and planet. It is given 
by 



L 



GMry 



-~~T COS (if - (f p ). 



(5) 



We adopt units such that G = M* = 1. The Keplerian or- 
bital period at r = 1 is then P(l) = 2tt. As in our previous 
work, for disc-planet simulations we use a uniform viscos- 
ity v — 10 -6 in code units. This is an order of magnitude 
smaller than the typically adopte d value of v — 10~ 5 that 
suppresses the vortex i nstability (|de Val-Borro et al ] l2007l : 
iLin fe Papaloizoull2010h . 



2.1 Disc-planet model 

Most of our discussion will be applied to disc-planet sys- 
tems, so we describe these models here. Supplementary sim- 
ulations in 3S1 employ a different set up from that described 
below. 

The disc occupies r = [n,r ] = [1,101- The initial 
surfac e density profile is modified from I Armitage &; Hansen! 
ET" 



E(r) = E r 



-3/2 



1 



n 



(6) 



where Hi = H(r\) is introduced to prevent a singular pres- 
sure force at the inner boundary because E — >• there. So 
is chosen via the parameter Q Q = QKe P (r ) where 



QKep(r) 



ttGE 7rr 2 E(r) 



(7) 



is the Toomre Q parameter for a thin Keplerian disc 
with a locally isothermal equation of state and k 2 = 
2Qr- 1 d(r 2 Q)/dr is the square of the epicycle frequency, with 
Q being the disc angular velocity. We use Q to label disc 
models and also define Q p = QKepfVp). We note that specify- 
ing Q also determines the disc-to-star mass ratio, M^/M*. 



The conversion between Q G , Q p and M^/M* for the mod- 
els used in disc-planet interactions is given by Table [ATI in 
Appendix [Al 

The disc is initialised with the azimuthal velocity re- 
quired by hydrostatic equilibrium. The initial radial veloc- 
ity is u r — 3v/r, corresponding to the initial radial velocity 
profile of a Keplerian disc with uniform kinematic viscosity 
and surface density oc r _3//2 . Note that \u r \ <^ \uip\- 

In most of this work we focus on stability of the gap in- 
duced by an embedded planet, accordingly we fix the planet 
on a circular orbit of radius r p = 5. We quote time in units 
Po = 2n /Qk(r p ). The planet is introduced at t = 25Po with 
azimuthal velocity that takes account the contribution from 
the gravitational force due to the disc. The planet potential 
is ramped up over a time period of 10Po- 



3 LINEAR MODES IN DISCS WITH 
STRUCTURED SURFACE DENSITY 

In this section and the next we study linear disturbances 
associated with internal surface density depressions in a 
disc in which self- gravity is not neglected. We will consider 
dips/gaps self-consistently opened by a giant planet. To sim- 
plify the discussion, we ignore viscosity, indirect potentials 
and the planet potential. The planet's role is then only to set 
up the basic state, assumed to be axisymmetric and defined 
by E(r), u^^r) and u r — 0. We begin with the linearised 
equations with a local isothermal equation of state. 

We consider perturbations to the disc state variables 
with azimuthal and time dependence through a factor 
expz(cr£ + mip) which is taken as read. Here a = cfr — ij 
is a complex frequency with <jr and —7 being the real and 
imaginary parts respectively. The azimuthal mode number 
m is taken to be a positive integer. Denoting perturbations 
by a prime, the linearised equations of motion read 



1_ 
D 



1 

D 



2 dW (W 

dr dr 

2 dW d& 

dr dr 
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2imQ f 2 

( c 
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ma 



c s W + $ x ) 



(c 2 a W + &) 



(8) 
(9) 



where D = k 2 — a 2 , a = a + raQ(r) is the Doppler shifted 
frequency, W = E x /E is the relative surface density pertur- 
bation and & is the disc gravitational potential perturbation 
which is given by the Poisson integral 



& = -G r° K m (r,Z)H(Z)W(t)£dZ, where (10) 
cos(mip)dip 



K, 



Jo 



^r 2 + £ 2 -2r£cos( V ) + e!(£)' 



(11) 



Using the linearised equations of motion to eliminate the 
velocity component perturbations from the linearised conti- 
nuity equation 



1 d , ,v im , 

icrW = — (rY>u r ) U4 

rE dr v J r * 
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dr 



+ 



+ 



yields the governing equation 
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a dr 
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W 



(c 2 s W + = 0. 



(13) 



3.1 Modes leading to vortex formation: the effect 
of self-gravity 

We shall be concerned with the vortex- forming instabili- 
ties commonly obse rved in disc-plane t interactions in discs 
with low viscosity jKoller et aljbooS iLi et al.1 12005L 120091 : 
iLin &; Papaloizoull2010h . The unstable modes are found to 
be localised to gap edges and associ a ted with minima in the 
vortensity (IPapaloizou Linl Il989l ; IPapaloizou &; Savoniid 
ll99ll : lLovelace et alJll999h . 

Here, we investigate the effect of self- gravity on these 
modes and show somewhat paradoxically, in view of the fact 
that increasing self- gravity in general destabilises discs, that 
that it tends to stabilise the modes. To do this we consider 
a small change to the strength of self-gravity and apply per- 
turbation theory to determine the consequence of the change 
for these modes. 

Non-linear simulations show more vortices develop as 
self- gravity is increased. This is in fact consistent with the 
stabilisation effect demonstrated below, which is effective 
for small m. In other words, stabilisation by self- gravity dis- 
courages low-ra vortex modes, supporting the observation 
that only higher m modes develop in simulations. 

To make the analysis tractable, further simplifications 
are made. We take the fluid to be such that either c 2 s is con- 
stant (strictly isothermal) or more generally the equation of 
state is barotropic such that p — p(S). We do not expect a 
significant difference between adopting a strictly isothermal, 
or barotropic equation of state, and the fixed c s profile as 
used in simulations. This is because the modes of interest 
are driven by local features and disturbances are localised, 
whereas c s (r) for a locally isothermal equation of state varies 
on a global scale. For mathematical convenience we also take 
a softening prescription such that Km(r,£) = K m (^r) is 
symmetric, eg. e g = constant. Although convenient mathe- 
matically, this is not expected to lead to significant changes 
for the same reason as that given above. 

Introducing the variable S 



S = c s W + 



c s W 



1 



G / JT m (r,0££(0W(0de. ( 14 ) 



the governing equation derived from equation ([T3|) simply 
by taking c s to be constant is 



dr 



77 

rL(S) 



dSX 
dr J 



+ 



m d 
a dr 



7,(1 



rD 



(15) 



where we have used the expression for vortensity 77 = 
k 2 /2QT i and v = a j k. We remark that in fact these equa- 
tions also hold for a general barotropic equation of state. 



3.2 The limit of negligible self-gravity 

When self- gravity is negligible, we may set G — in (|14[) 
from which we obtain W — S/c 2 s . Substituting this into ([T5]) 



gives the second order ordinary differential equation for S 
that governs stability in the limit of zero self-gravity in the 
form 



r^S 

r 2 



d_ 

dr 



D [dr 



^ m d 
a dr 



77(1 -v 2 ) 



m 2 S 
rD 



S. 
(16) 

This equation was studied by ILin &; Papaloizo"ul (|2010l ) for 
gap profiles of the type we consider. They found that de- 
pending on parameters such as m, there were neutral modes 
with co-rotation radius , r c , where a(r c ) = coincident with 
a vortensity minimum. This could be associated with either 
the inner or outer gap edge. Small variations of the disc pa- 
rameters could then lead to an instability associated with 
vortex formation. When unstable these modes retain r c to 
be close to the vortensity minimum with which they can be 
considered to be associated. 



3.3 Localised low m modes 

When m is small the neutral modes described above are 
localised around co-rotation and therefore insensitive to dis- 
tant boundary conditions. For these modes we may neglect 
v and set D - 



in ([16]) above to get the simpler equation 



r^S 

r-2 



dr 



rS (dSX 
dr~ ) 



+ 



m d 
a dr 



m 2 S 



S. 



(17) 



Localisation occurs because the solutions of (fTT)) can be seen 
to decay exponentially away from the co-rotation region, 
where there are large background gradients, on a length scale 
comparable to H. However, for (p"T[) to be a good approx- 
imation, we require |z> 2 | <C 1 in the region of localisation 
which is also comparable to H in extent. This in turn re- 
quires H <C 2r/(3ra), a condition which is satisfied for low 
m. For larger m the mode becomes de- localised with the ex- 
citat ion of density waves th at propagate into the extended 
disc ()Lin &; Papaloizoull2010h . In this case the boundary con- 
ditions can play a role. The analysis below assumes local- 
isation so that boundary conditions do not play a role. It 
therefore only applies for low m. 



3.4 Evaluating the effect of small changes to low 
m modes using perturbation theory 

In order to investigate the effect of self-gravity we return 
to equations (fT4)) and JTSJ) and note that strengthening self- 
gravity by scaling up the surface density is equivalent to 
increasing G, provided that the background form remains 
fixed. Thus although we increase the effective disc gravity, 
we shall assume the background disc model remains un- 
changed. In fact, in our case the background model is struc- 
tured by a perturbing planet and so its response to changing 
the disc gravity is non-trivial to evaluate analytically. Calcu- 
lations presented below in section 14.71 show that changes to 
the background surface density profile induced by incorpo- 
rating the disc gravity tend to act to make the vortex modes 
we consider more unstable. However, the direct effect of self- 
gravity on the linear response considered below turns out to 
be more important for localised modes with low m, and acts 
to stabilise them. 

Accordingly we assume a solution S corresponding to a 
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neutral mode with co-rotation radius, r c , located at a vorten- 
sity minimum, such that dr]/dr\ r=rc = 0. As the associated 
a is real, for this value, the operator L is real and regular 
everywhere. We then perturb this solution by altering the 
strength of self-gravity via 

G -> G + SG 

so that 5G > corresponds to increasing the importance 
of self- gravity and vice versa (note that the initial value, G, 
could be zero). This induces perturbations 

S + 5S, + a^a + Sa, 

L -> L + SL, 
with 



Sa : 
SL 



SCFR 

dL 
' da 



■ 27 



5a 



where Sa r and 7 are real. 

Noting that S denotes a small change and 7 is small, 
we linearise in terms of these quantities about the assumed 
original neutral mode and determine 7. The governing equa- 
tions lead to 



L(SS) + SL(S) = <SE' 



8S- 



G f r ° K m (r,t)£SE'(t)d£ 
SG f° K m (r,£)& '(Ode- 

J Ta 



(18) 



(19) 



Next, we define the inner product between two functions 
U(r), V[r) as 



(U,V) 



J Tn 



rU*(r)V(r)dr, 



(20) 



Then, assuming localised functions corresponding to lo- 
calised modes so that we can assume boundary values vanish 
when integrating by parts, we have 

(U,L(V)) = (L(U),V), (21) 

which used the fact that L corresponding to the neutral 
mode is real, which makes L self-adjoint. Equations II 8J— [T9l 
can be manipulated to yield 



(S,8L(S)) 



J Ta J Ta 



SGE, 



(22) 



where we have used the fact that the kernel K m is sym- 
metric. Note that E > and is proportional to the mag- 
nitude of the gravitati onal energy of the mode. Following 
IPapaloizou &; Linl (|l989h , we separate out the contribution 
to the perturbed linear operator SL that is proportional to 
the vortensity gradient and is potentially singular at co- 
rotation and write 



(S,5L(S)) = {S,5L!(S)) + (S,5L 2 (S)), 



(23) 



where 5L2(S) contains the potentially singular contribution 
and we have 



(S,SL 2 (S)) = -m8o 



r + r 

h % ° J It, 



where g(r) 



r dr 



m8(a)g(r)dr 
(24)' 

(25) 



Note that the usual Landau prescription for positive 7 has 
been used to deal with the pole at co-rotation and that 
V outside the above integral indicates that the principal 
value is to be taken. Note too that 5L± accounts for the 
remainder of SL and is not singular at co-rotation. The 
contribution from 8L2 is the only one that can lead to an 
imaginary contribution in the limit 5a —> because of the 
pole at co-rotation. In this regard we remark that one does 
not get such contributions arising from Lindblad resonances 
where D = 0. This is because, as is well known, these do 
not constitute effective sing ularities in a gaseous disc (eg. 
IPapaloizou Savoniiel[T99lh . 

Recalling that 7 is small, equations 1221- [24l can be com- 
bined to give 



SGE = Sa 



d&L^S)) 



da 

-im7r / S(a)g(r)dr 

J ri 

(A — imx)5a 



mV 



(26) 



where A is the contribution from the L\ term plus the prin- 
ciple value integral and 



X 



*\s\ 2 



mQ'lmQ'l dr 2 \rj , 
In the limit 7—^0+, A is real so we have 



(27) 
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A 2 + m 2 x 2 



SGE. 



(28) 



As we have remarked above, vortex forming modes are 
associated with vortensity minima or at maxima of 77 -1 . 
Consider a marginally stable mode with co-rotation at 
max(?7 _1 ). Typical rotation profiles have Q' < 0, which 
means x > f° r this mode. For consistency with the as- 
sumption of 7 > in deriving equation[28] we require SG < 
since E > 0. Accordingly in order to de-stabilise this mode, 
the strength of self- gravity needs to be reduced. 

This leads to the conclusion that: increasing self- gravity 
stabilises low ra modes with co-rotation at a vortensity min- 
imum (as has been borne out by our linear and nonlinear 
calculations presented below); while increasing self-gravity 
would de-stabilise modes which had co-rotation at a vorten- 
sity maximum. This suggests that for sufficiently strong self- 
gravity, modes associated with vortensity maxima should be 
favoured. This has been found to be the ca se but they are not 
mode s leading to vortex formation (see iLin &; Papaloizoul 
l201lh . 

Note x suggesting that 7 decreases for large 

ra, so the stabilisation/destabilisation effect of self- gravity 
diminishes for increasing azimuthal wave-number. This is 
only speculative because there are implicit dependencies on 
ra through terms in the integrals and through the original 
eigenfunctions S. Nevertheless, a weakening effect of self- 
gravity through the potential perturbation is anticipated be- 
cause increasing ra decreases the magnitude of the Poisson 
kernel K m - 
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3.5 Association of localised low m normal modes 
with vortensity minima for the strength of 
self-gravity below a threshold 

We can also show that corotation radii for localised neutral 
modes with low m must be at vortensity minima unless self- 
gravity ( or an appropriate mean value of Q) is above (be- 
low) a threshold level. Let us consider a disc with localised 
steep surface density gradients and a non axisymmetric dis- 
turbance with co-rotation radius in this region. Multiplying 
equation (|15[) by S* and integrating over the disc, assuming 
most of the contribution is from near corotation, as is ex- 
pected for low m modes ( see section I3.3[) , so that the term 
that is potentially singular at corotation and is proportional 
to the vortensity gradient is dominant on the right hand 
side, we have the balance 



3 E 

m\S\ 2 d (I 



rr 
J r,- 



~ 2 ' S ' -dr - G / K m (r,0r&*(0-Z'(r)dtdr 



a dr \T) 



- dr. 



(29) 



If the left hand side can be shown to always be positive 
then co-rotation for a localised neutral mode must lie at a 
vortensity minimum, or max(?7 _1 ). This holds if the maxi- 
mum possible value of A for any E' is less than unity, where 



A = 



G£°K m (r,t)r&*(£)l?(r)dtdr 



f^rc^dr 

Now from the Cauchy-Schwartz inequality it follows that 



(30) 



A < G 
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c2(r)c2(0 



r^d^dr. 
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Thus if the right hand side of the above is less than unity, 
corotation of a neutral mode localised at a vortensity ex- 
tremum must be localised at a vortensity minimum. This 
condition requires that the strength of self-gravity be below 
a threshold and this implies a lower bound on an appropri- 
ate mean Q value. The fact that this fails for sufficiently 
strong self-gravity is consistent with the discussion above 
that led to the conclusion that increasing self-gravity tends 
to stabilise the vortex forming instability. 

Indeed when self- gravity increases further, instability is 
transferred to modes with corotation associated with vorten- 
sity maxima. These are different in character to vortex form- 
ing modes being more global and are referred to as edge 
modes and they are studied in lLin &; Papaloizoul (|201lh . 



4 LINEAR CALCULATIONS 

We now present numerical solutions to equation (|13[) with a 
local isothermal equation of state for consistency with nu- 
merical simulations used to set up the basic state. We regard 
the governing equation as the requirement that an operator 
C acting on W should be zero, thus 



C(W) 0. 



(32) 



One form of C is 

f_ 

' dr 2 
'2m 
a 



C =rc s — — + rX 2 + 



1 T/_ _ U_ 
r + E ~ ~D 



rc s — + rXi 
dr 



+ 



2/ d 

+ rc s T + 
dr 



VLD' 
D 

-rD 



(33) 



It is understood that primes in equation l33l denote d/dr. 
The integro-differential operators X n are such that 



Xn(W) 



dr n 



(34) 



First and second derivatives of the perturbed potential are 
performed by replacing K m {r, £) by dKm/dr and d 2 Km/dr 2 
in the Poisson integral respectively. 



4.1 Quadratic approximation 

Vortex modes are localised about co-rotation at a vorten- 
sity minimum. They extend over a region characterised by 
small \a\ in the neighbourhood of r c and are expected to be 
insensitive to boundary conditions. If we multiply equation 
(|33[) by aD and expand the resulting equation in powers of 
a and neglect terms proportional to a 3 and higher powers, 
we can derive an equation of the form 



(cr 2 C2 + aCi + Co) W = 0, 



(35) 



where the operators d are real and only depend on the 
basic state. The eigenfrequency now appears explicitly. We 
call this the quadratic approximation. Although the above 
procedure is not expected to be valid in general, it should be 
valid for modes localised about their corotation radii with 
the scale of localisation being much less than r c itself. Modes 
can be found from ([35]) using standard numerical methods. 
The eigenvalues are then used as starting values in order to 
obtain an iterative solution for the eigenvalues for the full 
equation ( [32]) . The concept of the vortex forming mode as 
a localised mode is confirmed if the final solution of ( [32]) is 
not significantly changed from the solution of (|35[) . Note that 
since the operators & are real, eigenfrequencies are either 
real or occur in complex conjugate pairs. 



4.2 Distinction between including and excluding 
self-gravity 

Since we are concerned with effects of self-gravity, We 
need to clearly define self- gravitating (SG) and non-self- 
gravitating cases (NSG). The background state on which we 
perform linear stability analysis were obtained from non- 
linear hydrodynamic simulations (see ^4.4|) . which can be 
run with self-gravity (SGBG) or without self-gravity (NS- 
GBG). In the linear calculations, self-gravity can be dis- 
abled by setting $ ; = 0. Therefore we have two disturbance 
types: a self- gravitating response (SGRSP) and a non-self- 
gravitating response (NSGRSP). 

The analytical discussion in 33.11 applies to the effect 
of self- gravity through the linear response, assuming the 
form of the background remains unchanged. In the con- 
text of hydrodynamic simulations, including or neglecting 
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self-gravity simply means whether or not the disc gravi- 
tational potential is included. If it is included, the back- 
ground state set up by simulations also depends on self- 
gravity. Hydrodynamic simulations therefore correspond ei- 
ther to the combination SGBG+SGRSP or to the combina- 
tion NSGBG+NSGRSP. We call these fully self-gravitating 
and fully non self- gravitating cases. The advantage of per- 
forming linear calculations is that we can distinguish be- 
tween the effects of self-gravity arising through its effects 
on the background state and the effects resulting from its 
influence on the linear response. 



4.3 Numerical approach 

We discretised the linear operators on a grid that divides the 
radial range [ri,r ] into typically N r = 385 equally spaced 
grid points at which W is evaluated as Wj, j = 1, 2..N r . The 
governing equation thus becomes a system of N r simulta- 
neous equations for the Wj which determine an eigenvalue 
problem for a. The system of equations incorporates the 
boundary conditions. For simplicity we impose dW/dr — 
at boundaries in discretised form. We remark that it is 
known from other linear calculations and simulations that 
vortex-forming modes are localised and insensitive to bound- 
ary c onditions (|de Val-Borro et all 120071 : iLin &; Papaloizoul 
2010). Tests have shown that the boundary condition does 
not influence the essential effect of self-gravity on localised 
modes. 

The discretised quadratic approximation obtained from 
equation (|35[) leads to a quadratic eigenvalue problem that is 
equivalent to a 2N r x 2N r standard linear matrix eigenvalue 
problem. This is solved by methods which give all eigenval- 
ues a. The most unstable eigenvalue is then used as a trial 
solution for the full system which yields a discriminant that 
is then solved iteratively for the actual eigenvalue , a, us- 
ing the Newton Raphson method. We restrict attention to 
modes with eigenfrequencies such that \<jR/mQ e + 1| < 0.1 
and I7I < Q e , where Q e is the rotation frequency of the 
vortensity minimum at the outer gap edge. 



4.4 Background state 

The background states used for linear calculations were set 
up by running disc-planet simulations for the models de- 
scribed in 32.11 Numerical details will be described in JSl 
We allow the planet to open a gap, then take azimuthal 
averages to obtain one-dimensional profiles. 

Fig|T] shows the gap profile, in terms of the inverse 
vortensity I/77, opened up by a Saturn-mass planet in the 
Qo = 4 disc. The formation of vortensity peaks is due to the 
generation of vortensity as material passes through shocks 
induced by the giant planet whi ch extend into its co-orbital 
region (|Lin fe Papaloizoul [2010). We shall see that in hy- 
drodynamic simulations vortices predominantly develop at 
the outer gap edge, so we focus on modes associated with 
the outer vortensity extrema. The outer vortensity maxi- 
mum or min(?7 _1 ) is located at r ~ 5.5 while the outer 
vortensity minimum or max(?7 _1 ) is located at r c± 5.75. 
These extrema are separated by about 0.89if . As we in- 
crease the strength of self-gravity, the discussion in sections 
13.41 and 13.51 indicates that eventually vortex modes associ- 
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Figure 1. Gap profiles opened by a Saturn-mass planet in an 
initial disc model with Q Q = 4. The inverse vortensity r]~ 1 , ob- 
tained with (solid) and without (dotted) self-gravity included is 
shown. Profiles were obtained from azimuthal averages taken from 
disc-planet simulation outputs. The planet is at the fixed location 
r = 5. 

ated with vortensity minima are suppressed and modes as- 
sociated with vortensity maxima become favoured instead. 
This is consistent with self-gravity stabilising low m vortex 
modes through the linear response. 

The gap profiles set up with and without self-gravity 
are similar. When self- gravity is included, the peaks and 
troughs have slightly larger amplitudes due to increased ef- 
fective planet mass (see fc|5.3[) . 

4.5 Solution in the quadratic approximation 

We first present solutions for linear modes for the Q = 4 
disc in the quadratic approximation. Fig. [2] compares the 
growth rates ,7, and the eigenfunctions, W, obtained for 
m = 5 for the fully self- gravitating case and the fully non 
self- gravitating case. Growth rates are such that the most 
unstable mode shifts to higher m when self-gravity is in- 
cluded. Without self- gravity, the dominant mode has m = 3, 
whereas with self-gravity, the dominant mode has m — 6 — 
7. The combination of self-gravity acting through the back- 
ground and the linear response stabilises modes with m ^ 5 
and de- stabilises modes with m ^ 6. Thus higher m vortex 
modes are enabled by self-gravity. 

The eigenfunctions W with and without self-gravity 
are similar, and are for the most part localised about co- 
rotation. This suggests the nature of the instability is un- 
changed in this case. However, the eigenfunction amplitude 
inr>6,r<5.2is larger when self- gravity is included. 
Note too that the amplitude is larger for r > 6 than for 
r < 5.2. This is not unexpected since the outer disc has 
smaller Toomre Q values it is expected to be more respon- 
sive. 

It is important to remember that the quadratic approx- 
imation assumes modes are localised about co-rotation. The 
global nature of some disturbances may invalidate this ap- 
proximation. Increasing m eventually quenches modes in the 
non self- gravitating case, thus we did not find modes for 
m ^ 9 that fit the description of being localised modes. For 
the full governing equation (fT3)) , high m modes are expected 
to have increasing wave-like behaviour and be less focused 
around corotation, contradicting the quadratic approxima- 
tion and so it is not surprising that they are not found here. 



8 Lin and Papaloizou 



2.5 














(a) ; 


2.0 








1 .5 








1 .0 




















NSGBG+NSGRSP 




0.0 









1 2 3 45 6 7 8910 



m 

rn= 5 



1 .2 






A (b) - 


1 .0 






0.8 




NSGBG+NSGRSP 




0.6 








0.4 








0.2 








0.0 









2 3 4 5 6 7 8 

r 

Figure 2. Growth rates (a) obtained from the quadratic approx- 
imation to the governing equation and the modulus of the m = 5 
eigenfunction (b). Here W c = W(r = 5.7). The solid (dotted) 
lines correspond to the case with (without) self-gravity in both 
the background state and response. The initial disc model had 
Qo = 4. 

Similarly, we did not find a localised m — 1 vortex mode 
in the self- gravitating case because it had been stabilised 
by the inclusion of self- gravity, according to earlier analysis. 
However, there exists other types of low m mode not cap- 
tured by the quadratic approximation. Modes with extreme 
values of m are of less relevance since they are not observed 
to develop in hydrodynamic simulations of interest, which 
typically show the number of vortices in agreement with the 
most unstable m. 



4.6 Solutions to the full governing linear equation 

Solutions to the full governing equation (fT3|) , for the fidu- 
cial cases with Q Q = 4 above are shown in Fig. [3] In Fig. [3] 
we compare growth rates and the m = 5 eigenfunction. The 
corotation radii for the self- gravitating and non-self gravi- 
tating cases are r c = 5.88, 5.79, respectively. These radii are 
close to the local vortensity minimum or max(?7 _1 ). In the 
non self- gravitating case, co-rotation almost coincides with 
the extremum. 

The non- self- gravitating case has maximum growth rate 
for m — 4 and the self- gravitating case is most unstable for 
m = 6 — 7. Modes with m < 5 are stabilised while m > 5 
are de-stabilised by self-gravity. This results in a shift to 
higher m modes when self-gravity is included in nonlinear 
simulations of the model. Modes with the larger m values are 
destabilised, but this is not in contradiction with with earlier 
analysis which assumed small m and did not take account of 
variations in the form of the background state. It is shown in 



section [47T1 below that the destabilisation of higher m modes 
is in fact due to the change in the background state. 

Results for m < 5 are in essential agreement with those 
obtained in the quadratic approximation giving us confi- 
dence in our view of the importance of the dominance of 
the corotation region. As with the quadratic approximation, 
we did not find a m = 1 mode with dominant disturbance 
about the gap edge in the self- gravitating case. This suggest 
such modes are the most easily stabilised by increasing self- 
gravity. In the non self- gravitating case, the growth rates for 
m ^ 8 are larger than m — 1 and do not follow the trend of 
decreasing growth rates seen from m — 5 —> 7. Notice in Fig. 
[3ja) (and Fig. Ufa) below) the local maxima at m = 9 for 
non- self- gravitating cases. The non-self-gravitating m ^ 8 
modes are unlikely to be the same type of vortex modes as 
m ^ 7 because the former have significant wave-like regions 
in W and are not concentrated near co-rotation as for m ^ 7. 
In addition, the m ^ 8 region is also where the quadratic ap- 
proximation appears to fail. These wave- dominated modes 
demand radiative boundary conditions rather than the sim- 
plistic conditions applied here, which are appropriate for 
boundary condition insensitive modes such as local vortex 
forming modes. The non-self- gravitating m ^ 8 modes iden- 
tified here are thus likely to be artifacts of the boundary 
condition. Fortunately, these modes are irrelevant because 
they are not the most unstable, nor are they observed in the 
corresponding nonlinear simulations. By considering the be- 
haviour for m ^ 7, we can expect a cut-off for vortex modes 
around m = 8. for this model. 

The m — 5 eigenfunctions in shown in Fig. [3] are sim- 
ilar. Both have dominant disturbance around co-rotation. 
Behaviour in the region r ^ 4.6 is essentially identical. The 
largest difference is found in the region r ^ 6.4. The dis- 
turbance around co-rotation is also somewhat wider in the 
self- gravitating case. Comparing with Fig. [2j we see that 
the quadratic approximation captures the main feature of 
the mode in the co-rotation region. However, it removes the 
wave-like behaviour in the exterior disc. 

4.7 The role of self-gravity 

The calculations above can be compared to hydrodynamic 
simulations where self- gravity is either enabled or not. The 
linear problem allows one to examine separately the effect 
of self-gravity through its influence on the basic state and 
through its influence on the linear response. 

We continue with the disc model with Q Q = 4 and com- 
pare growth rates for a pair of models, one with and the 
other without self-gravity in setting up the basic state (i.e. 
the simulations were run with and without self-gravity en- 
abled) but both without self-gravity implemented in the lin- 
ear mode calculation (i.e. setting 4>' = 0). In addition we 
compared a another pair of models, one with and the other 
without self-gravity implemented in the linear response, but 
with both having the background state calculated with the 
disc self- gravity incorporated. 

Fig. [HJa) examines the influence of self-gravity through 
its modification of the background state. It shows that this 
modification is de-stabilising. This is not surprising because 
a deeper gap is set up when self-gravity is included in the 
simulation. In going from the case without self- gravity to 
the one with self-gravity, the most unstable mode shifts from 
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Figure 3. Growth rates obtained for the full governing equation 
(a) and the eigenfunction \W\ for m = 5 (b). Here W c = W(r = 
5.7). Solid lines are for the fully self- gravitating case and dotted 
lines are for the fully non self-gravitating case. The initial disc 
had Q = 4. In (a), the local maximum around m = 9 in the 
dotted curve is most likely caused by a boundary condition effect. 
Modes with m ^ 7 without self-gravity are not seen in nonlinear 
simulations and are not relevant. 



m — 4 — 5 to m = 6. The vortex modes do not require self- 
gravity to operate. However, the difference in growth rates 
decreases as m decreases so the effect due to the modification 
of the background is smallest for small m. In other words, 
high m modes can be made more unstable by including self- 
gravity in the base state. 

FigHJb) compares the growth rates obtained from lin- 
ear calculations for the same background state but with and 
without self- gravity implemented in the response. The anal- 
ysis in A3. II applies to this comparison. Consistent with that, 
enabling self- gravity in the response decreases \~f\ and sta- 
bilises the system against modes with co-rotation associated 
with a vortensity minimum for any value of m. Unlike the 
effect acting through the background, the influence of self 
gravity is more significant for low m and can lead to sta- 
bilisation. The m — 1 mode has been stabilised. The most 
unstable mode without self-gravity has m = 5 and including 
self- gravity shifts it to m = 6. 

The fact that self-gravity acting in the linear response 
and background state has opposite effects on growth rates is 
consistent with the comparison between fully self- gravitating 
and fully non-self-gravitating cases (Fig. [3}. The effect of 
self-gravity through changes to the background and through 
direct influence on the linear response both contribute 
to favouring higher m. However, the background effect is 
achieved by increasing high m growth rates, whereas the ef- 
fect via the response works by stabilising low m modes in 
accordance with the discussion in section 13.41 Thus low m 



Figure 4. The effect of self- gravity on growth rates (as a function 
of azimuthal wavenumber m) of vortex-forming, gap edge instabil- 
ities through the background (a) and through the response (b). 
The disc model is Q = 4. In (a), the local maximum around 
m = 9 in the dotted curve is most likely caused by a boundary 
condition effect. Modes with m ^ 7 without self- gravity are not 
seen in nonlinear simulations and are not relevant. 



vortex modes become disfavoured. Overall, we expect more 
vortices to form, corresponding to increasing m, with in- 
creasing self-gravity. 



4.8 Models with different Q Q 

In the calculations presented below, self-gravity is included 
both in setting up the background and in evaluating the 
linear response. We compare the fiducial calculation to cases 
with Q = 3 (a more massive disc) and Q — 8 (a less 
massive disc). The Q Q = 8 case has a mass Md = 0.012M* 
which is usually considered to be non- self- gravitating in disc- 
planet simulations. 

Growth rates are shown in Fig. [5{a). The plot demon- 
strates a clear shift of the most rapidly growing mode to 
higher m as Q is lowered (increasing disc mass). For Q Q — 8 
the most unstable mode has m = 5 — 6 and for Q Q = 3 it 
shifts to m = 7 — 8. The shift is accompanied by the stabili- 
sation (or loss) of low m modes. For example, when the disc 
mass is halved as Q Q changes from Q Q = 4 to Q Q = 8, the 
m — 3 growth rate decreases by a factor of 3. Modes with 
the lowest m are stabilised by strong self-gravity, as we did 
not find m = 1, 2 modes when Q Q = 3. 

These calculations have also been performed in the 
quadratic approximation. Fig. [5jb) shows a similar depen- 
dence of the growth rates on Q Q in this case. The shift to 
higher m is more apparent. For Q Q = 3 we did not find 
modes with m ^ 4. Since the approximation reinforces the 
localised property of vortex- forming modes, it means that 
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Figure 5. Growth rates for unstable modes in background discs 
with Q = 3, 4, 8. (a) is obtained from solving the full linear 
equation, while (b) is from the quadratic approximation. Self- 
gravity is fully incorporated throughout. 



localised modes for low m become increasing unlikely as we 
lower Q . Hence, for massive discs only high m vortex form- 
ing modes can develop. The agreement between results ob- 
tained using the quadratic approximation and the full gov- 
erning equation indicates a lack of sensitivity to boundary 
conditions and so is reassuring. 



5 NONLINEAR HYDRODYNAMIC 

SIMULATIONS OF VORTEX FORMING 
INSTABILITIES 

We present hydrodynamic simulations of vortex formation 
and evolution at edges of gaps opened by a giant planet 
for discs of varying masses with self-gravity self-consistently 
included. The disc-planet models have been already de- 
scribed in £12.11 We consider a Saturn-mass planet with 
M p = 3 x 10~ 4 M*. The kinematic viscosity v = 1CT 6 cor- 
responds to an a viscosity of a = v/(c s H) — 1.8 x 10 -4 at 
the planet's fixed orbital radius at r p — 5. 

Higher viscosities corresponding to a = O(10 -3 ), which 
is commonly as sumed for protoplanetary discs, inhibit vor- 
tex f ormation (|de Val-Borro et all 1 20071 ; iLin &; Papaloizou! 
l201Ch . If planetary masses appropriate for type I migration 
are used, then a much lower visco sity would be required 
for vortex formation to occur (e.g. iLi et all 12009. where 
a ^ 10 -5 was needed). Since we consider a giant planet, 
such low viscosity regimes are not needed in order to de- 
velop vortices. 



5.1 Numerical scheme 

The h ydrodynamic equations are evolved using the FARGO 
code (Masset 2000). F ARGO is an explicit fi nite- difference 
code similar to ZEUS (|Stone fe Norman| [l992) with second 
order accuracy in space. It employs a modified azimuthal 
transport to achieved large time-steps, otherwise limited by 
the Keplerian velocity near the inner boundary of the do- 
main . Self- gravity for FARGO was implemented and tested 
by iBaruteau &; Massetl (|2008[ ) . The gravitational accelera- 
tion due to the disc is calculated directly using Fast Fourier 
Transforms in both azimuth and radius. The latter requires 
the radial domain to be doubled when calculating the self- 
gravity potential. 

The disc is divided into NrXN^ = 768 x 2304 computa- 
tional grid cells in radius and azimuth. The computational 
grid in the radial direction is logarithmically spaced. We im- 
pose an open boundary at r = n and non-re flecting bound - 
ary as used bv lZhang et~ai1 <|2008h at r = r Q (lGodonlll996h . 
Since vortices are localised features, as long as gap edges 
are far from boundaries, boundary conditions can only have 
a limited effect. We also perfor med some simulations with 
damping boundary conditions (|de Val-Borro et~al ] l2007h 
and open outer boundaries. We found similar results to those 
presented below. 

5.2 The effect of self-gravity 

We compare cases where disc gravitational potential is ei- 
ther included or excluded in the total potential calculation. 
This is the standard distinction between self-g ravitating and 
non-self-gravitating disc- planet simulations ([Nelson fc Bend 
l2003al : IZhang et al.ll2008h . We consider the Q Q = 4 ( Q p = 7) 
disc. The disc mass is M d = 0.024M*. 

Fig. [6] shows vortensity contours that illustrate vortex 
formation at the outer gap edge. There is a local vortensity 
maximum at about r = 5.5 that has been produced by flow 
through shocks induced by the planet. This maximum re- 
mains largely undisturbed indicating that instability is as- 
sociated with structure outside it and associated with the 
vortensity minimum. More vortices are excited when self- 
gravity is included. In that case the m — 6 vortex mode 
dominate whereas the m — 3 mode dominates in the non 
self- gravitating case. This behaviour is consistent with lin- 
ear calculations. 

Fig. [6] also shows that vortices have radial length-scales 
comparable to the local scale- height (~ H(6) = 0.3). The 
vortices in these cases have radial sizes of ~ 2.2if (6) without 
self-gravity and ~ 1.3if (6) with self-gravity. The vortices are 
of smaller radial extent when self-gravity is included because 
of the preference for higher m. A decrease in radial size 
can be expected from the decrease in width of the WKBJ 
evanescent zone centred on corotation. This is determined 
by the condition 

O + mQ) 2 = k 2 (1 - 1/g 2 ). (36) 

Writing a — -mQ(r c ), it is straight forward to show that the 
radial width of the evanescent zone approximately scales as 
1/m. Since the preferred m is approximately a factor of two 
smaller, we expect vortices in the case without self- gravity 
to be double the size of those in the case where self-gravity 
is included. 
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Figure 6. Vortensity contours for Q = 4 with (left) and with- 
out (right) self-gravity. The planet is located at r = 5, <j> — 4>p- 
The thick lines crossing the outer radial boundaries correspond 
to spiral shocks induced by the planet. 

In the case with self-gravity, the vortices are approxi- 
mately centred along the radius r = 5.9, close to the co- 
rotation radius expected from linear calculations r c = 5.88. 
In the non self- gravitating case, linear calculation gives 
r c = 5.81 for m = 3, but the vortices are approximately cen- 
tred along r = 6. However, perturbations here are already 
in the non-linear regime and interaction between vortices or 
with the spiral shock may shift the vortices around. Note 
too that in this regard there is more variation in the ra- 
dial locations of the vortices as compared to the case with 
self- gravity. 

5.3 Varying disc mass: gap profiles 

We present simulations of discs models with 1.5 ^ Q ^ 8, 
equivalents 2.6 < Q p < 14 or 0.063 ^ M d /M* ^ 0.012. The 
equilibrium gap profiles have a range of Q values with local 
extrema of 1.5 < min(Q) ^ 9.5 and 4.8 < max(Q) < 21.6 
near the outer gap edge. The gap profiles opened by the 
planet are given in Fig. for a range of Q that develop 
vortices. Outside the plotted region the curves are indis- 
tinguishable. The gap deepens with decreasing Q Q and gap 
edges become steeper. In going from Q Q = 8 —> 2 the gap 
depth |AE/E| increases by about 0.05 — 0.08, similarly the 
bumps near gap edges increase by 0.05 — 0.06. On a global 
scale, self- gravity becomes more important with increasing 
radius. However, we observe no trend in the difference be- 
tween gap profiles with respect to radius. 

Self-gravity affects gap structure on a local scale by in- 
creasing the effective planet mass so that M p — > M' p . A 
straightforward estimate, based on the unperturbed disc 
model, of the expected mass within the Hill radius , rjj, of a 
point mass planet with M p = 3 x 10~ 4 M* is M H 0.047M P 
for Q = 8 and M H = 0.17M P for Q = 2. It is likely that 
Mh , or at least some significant fraction of it adds to the ef- 
fective mass of the planet acting on the disc when self-gravity 
is self-consistently included. Thus M p , and therefore also 
the gap depth, is expected to increase with disc mass. Thus 
we note that without carefully tuning M p , self- gravitating 
disc-planet calculations with different surface density scales 





0.2 










0.1 






It' '^^fciw - 




0.0 






o 


-0.1 




Go =2 \ 






-0.2 




o = 3 \ //■ /- 

o =4 Y/ 

Q o =8 ^ 






-0.3 









3 4 5 6 7 



Figure 7. Gap profiles opened by a Saturn mass planet in self- 
gravitating discs as a function of disc mass, parametrised by Q - 
The azimuthally averaged relative surface density perturbation is 
shown. The planet located at r = 5. 

will always differ in M' p . Since the gap profiles differ, as 
we have seen above, stability is affected also. Note that the 
gap width w in Fig. [7] does not change greatly with Q , 
which is consistent with the scaling w oc m oc M p 1/3 . How- 
ever, the peaks/troughs for Q Q = 2 do lie slightly closer 
to the orbital radius r p than other cases. This is because 
shocks are induced close r to the planet due to increased M' p 
(jLin fe Papaloizoull20Ioh . 



5.4 Varying disc mass: gap stability 

Fig. [8] shows snapshots of the relative surface density per- 
turbation as instability sets in. Consider Q Q ^ 2 first. The 
instability is associated with the outer gap edge while the 
inner edge remains relatively stable. In the least massive disc 
Qo = 8, 3 — 4 vortices form at the outer gap edge, similar to 
what happens in the non- self- gravitating disc in £|5.2I As we 
increase the disc mass, more vortices develop. By Q Q — 2, 8 
surface density maxima can be identified. Note that a vortex 
may be obscured if it coincides with the planetary wake. In 
moving from Q Q = 8 — >> Q Q = 2, vortices become smaller 
and their centres move inwards. When Q Q = 8 local surface 
density maxima lie just outside r = 6 while for Q Q = 2 they 
lie just interior to r = 6. 

Increasing Md makes the perturbation more global. Vor- 
tices in massive discs can gravitationally perturb parts of the 
disc further out, similar to a planet. Lowering Q G , vortices 
develop longer, more prominent trailing spirals exterior to 
them. This is most notable with Q Q = 2, where the vortex 
spirals can have comparable amplitudes as the planet wake. 
However, increasing self- gravity even further, we expect a 
global instability to eventually develop as the vortices per- 
turb the disc strongly via self-gravity. 

The Qo = 1.5 disc does not develop vortices. This 
is consistent with the stabilisation effect of self-gravity on 
the vortex instability through the linear response. Instead, 
Qo = 1.5 develops global spiral instabilities associated with 
the gap edge. We call these edge modes. They have been 
suggested to oc cur in self-gravitating discs w ith surface den- 
sity depression (|Meschiari &; Laughlidl2008h . Here, we have 
shown that they can indeed develop in gaps self-consistently 
opened by a giant planet. Detailed discussions of edge modes 
are beyond the scope of this paper and are presented in 
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Figure 9. Fourier amplitudes of the surface density in the outer 
disc r E [5, 10], normalised by the m = component, for some of 
the models illustrated in Fig. [8] 



iLin &; Papaloizoul (|201lh , but for reference we note some im- 
portant differences to the vortex modes. Edge modes here 
are associated with local vortensity maxima, which lie closer 
to the planet than vortensity minima. Perturbations at the 
inner gap edge are also identified and correlate with those 
at the outer gap edge, where the over-density is largest (and 
stronger than for vortex disturbances). Communication be- 
tween the two sides of the gap is only possible with sufficient 
self- gravity. 

The Fourier amplitudes of the surface density in r G 
[5, 10], as a function of m, is shown in Fig. [9] for Q Q = 2, 4, 8 
at t = 56Po- Amplitudes have been normalised by the ax- 
isymmetric component. The shift to higher m vortex modes 
with increasing disc mass is evident as expected from lin- 
ear calculations. For Q Q = 8, 4, 2 the preferred vortex modes 
have respectively m — 4, 5 — 7 and 7 — 9, with average am- 
plitudes that decrease with Q . The latter may reflect the 
stabilisation effect of self-gravity on linear modes with low 
m. 

An important region in Fig. [9] is m ^ 3. There is a peak 
in amplitude at m — 2 for Q Q = 2, 4. These are the edge 
modes described above. They were not seen in linear calcu- 
lations because there we focused on finding vortex modes, 
which have co-rotation at or close to local vortensity min- 
ima. The loss of low m vortex modes with increasing self- 
gravity observed in linear calculations, is replaced by the in- 
creasing prominence of global edge modes. Fig. [9] shows the 
m = 2 amplitude becomes more significant with increasing 
self- gravity. 

In the Q = 2 simulations we do see evidence of an 
m — 2 disturbance hindering vortex evolution. This transi- 
tion case is difficult to analyse since both types of instabil- 
ities develop. Fig. [9] shows the edge and vortex modes have 
comparable amplitudes in r > 5. However, considering the 
region r ^ 7 for Q Q = 2 we found m — 2 is dominant, 
because the edge mode is global whereas vortex mode dis- 
turbances are localised to the edge (r c± 6). Increasing self- 
gravity further, we expect eventually edge modes become 
dominant, this is seen in the Q Q — 1.5 case in Fig. [8] 



5.5 Evolution and merging of self-gravitating 
vortices 

We examine the evolution of vortices under the influence of 
self-gravity. Fig. [10] compares surface density perturbations 
for a range of disc masses at t — lOOPo. Since the vortices 
emerge at roughly t = 56Po , they have evolved for a similar 
time. Typical growth rates for vortex modes are r ~ 6Po, 
which implies the vortices have evolved for about 7r, well 
into the non-linear regime. 

The snapshot for Q Q = 8 shows a single vortex, resulting 
from the merging of the initial vortices. The vortex distur- 
bance is largely confined to within a local scale-height of the 
gap edge . Such a result is typical fo r simulations with no self- 
gravity (|de Val-Borro et al]|2007l V The case with Q Q = 4 
shows vortex merging taking place, as individual surface 
density maxima can still be identified in the large vortex 
behind the shock. Notice there is a smaller vortex just pass- 
ing through the shock. When Q Q = 4, a single vortex forms 
at t ~ llOPo. Fig. [10] shows that increasing the disc mass 
delays vortex merging. For Q Q = 3.5, 5 vortices remain and 
for Q = 3 and Q = 2.5, 7 vortices remain. They have not 
merged into a single vortex as happened when Q Q = 8. A 
25% increase in disc mass as Q = 4 — >> 3 causes merging 
to be delayed by 50Po-This suggests that increased gravita- 
tional interaction between vortices opposes merging. 

As Q is lowered the inter-vortex distance increases. 
When Q — 2.5 their azimuthal separation can be larger 
than the vortex itself. Also, vortices become less elongated 
and more localised being symptomatic of gravitational con- 
densation. Trailing wakes from vortices also become more 
prominent as the vortices more strongly perturb the disc via 
their self- gravity. In fact, the planetary wake becomes less 
identifiable among the vortex wakes. Notice vortex wakes 
are mostly identified with vortices upstream of the plane- 
tary wake, rather than just downstream. The vortex wakes 
appear to detach from the vortex after passing through the 
planetary shock. There is a sharp contrast in gap structure 
between the Q — 2.5 and Q = 8 cases. The single vortex 
that results when Q Q = 8 is aligned along the outer gap 
edge, which is still approximately identified as circle r = 6. 
However, when Q Q — 2.5 the vortices and their wakes inter- 
sect the circle r = 6 making the radius of the outer gap edge 
less well-defined. 

An important feature of vortices is that the wakes 
are associated with the excitation of density waves which 
transport angular momen tum transport outwards (see e.g. 
IPaardekooper et "all l2010h . A measure of this is is made 
through the standard a viscosity parameter. We define it 
here through 



a(r,ip) = Au r Au(p/c^(r) (37) 

where A denote deviation from azimuthal averaged values. 
As we are interested in the vortex region, a is spatially av- 
eraged over the annulus r £ [5.7, 7.1] and its running-time 
average (a) is plotted in Fig. 1111 The parameter (a) associ- 
ated with vortices is O(10 -3 ), an order of magnitude larger 
than the imposed value associated with v. When Q Q = 4, 
(a) decays steadily after an initial transient growth. The 
case Q = 2 is also shown. In that case (a) ^ 1.9 x 10 -3 for 
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Figure 8. Instability at the outer gap edge of a Saturn-mass planet, in a discs with minimum Toomre parameter, from left to right, of 
Qo = 1.5, 2, 3, 4, 8. The relative surface density perturbation is shown. 
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Figure 10. Non-linear evolution of vortex instabilities at the outer gap edge of a Saturn mass planet as a function of disc mass, 
parametrised by the minimum Toomre parameter Q Q , from left to right, of Q Q = 2.5, 3.0, 3.5, 4.0, 8.0. 



t > 80Po- Interestingly for Q Q = 2.5 — 3.5 there is growth in 
(a) over several tens of orbits. 

The behaviour of (a) as a function of Q is consistent 
with the general picture of vortex formation and evolution 
(Fig. I10p . We see below that the decay in (a) is associated 
with vortex merging leading to fewer vortices. This is the 
case for Q Q = 4. Vortex merging happens more readily for 
lower disc masses, hence although multiple vortices develop 
from the instability, this phase does not last long enough for 
the multiple-vortex configuration to significantly transport 
angular momentum. As we increase self-gravity when Q Q = 
3.5 — ► 2.5, vortices become less prone to merging. Hence, 
the multiple- vortex phase lasts longer. They have time to 
evolve into compact objects that further perturb the disc. 
This is consistent with fact that (a) -growth is prolonged 
with increased disc mass, as merging is delayed. However, if 
self-gravity is too strong, such as when Q Q = 2.0, growth is 
again limited, because the m = 2 edge mode develops and 
hinders vortex evolution. 



5.6 Long term evolution of gap edge vortices 

We consider the long term evolution of vortices in the disc 
with Qo = 3. Fig. [12] shows the instantaneous a viscosity 
parameter, average over- density (defined in regions where 
the relative surface density perturbation is positive) and 
surface density contour plots of the vortex region. The vis- 
cosity parameter a grows from t = 100Po to t = 170Po 
with max(a) ~ 8 x 10 -3 , which is 50 times larger than the 
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Figure 11. The running time average of the viscosity parameter 
a averaged over the vortex region for a range of disc models. 



contribution from the background viscosity. At t = 150Po 
there remains 6 distinct vortices, one of which just passing 
through the planetary wake. The multiple- vortex configura- 
tion has been maintained for a further^ 50Po since Fig. 1101 
With weak self-gravity, a single vortex would have formed 
through merging. Delayed merging allows individual vortices 
to evolve and become planet-like. The typical over-density 
in a vortex is ~ 1 at t = 100Po and increases to > 1.5 by 
150Po- As a consequence we expect increased disturbance in 
the surrounding disc. The phase of a-growth correlates with 
a linear increase in the average over-density in the vortex 
region. 

At t = 170Po, Fig. n~2l shows 6 vortices still remain, 
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Figure 13. Comparison of post-merging vortices for discs with 
Qo = 3 (left) and Q Q = 4 (right). The Toomre Q (solid) and 
surface density (dotted) averaged over r G [6,6.5] is shown as a 
function of azimuth relative to the planet. 
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but the pair just upstream of the planet is merging. The 
parameter a is a maximum. However, after a burst of vortex 
merging events, a decreases rapidly. The snapshot at t = 
I8OP0 shows a quieter disc with only 3 vortices of similar 
size to those before merging. This is unlike cases with weak 
self-gravity where a larger vortex results from merging. Fig. 
1131 compares post-merging vortices in discs with Q Q = 3 and 
Qo = 4. For Qo = 3, the post-merging vortices are localised 
in azimuth, with Q ~ 1 and the contour plot shows the 
densest vortex has an over- density of ~ 3.75. However, for 
the disc with Q Q = 4, a single vortex forms that extends 
about half the total azimuth and has Q > 2. . 

The mass of the disc with Q Q = 3 is Md — 0.031M*, 
usually considered insufficient for gravitational collapse by 
itself. However, the comparison above shows that addition 
of self-gravity to the vortex instability induced by an embed- 
ded planet, may enable collapse into compact objects that 
survive against shear. 

We found that vortices in the disc with Q Q — 3 notice- 
ably affect the gap structure. Fig.ll4lshows several snapshots 
of the gap structure from t = 100Po, when there were multi- 
ple vortices, to the end of the simulation when a vortex pair 
remained. The gap profile for Q Q = 4 is also shown for com- 
parison. Neither the single vortex for the disc with Q Q — 4 
nor multiple vortices with Q Q — 3 affect the one-dimensional 
gap profile at t — 100Po because the disturbances only re- 
distribute mass in the azimuthal direction. However, the 
original bump at the outer edge is diminished after vortex 
merging takes place in the disc with Q Q — 3 at t — 200Po. A 
surface density depression of AE/H ~ —0.1 then develops at 
r = 7.3 and a bump develops at r = 8.5. These features last 
to the end of the simulation. Self- gravitating vortices behave 
like planets and gap opening is to be expected. Assuming 
vortices lie near the surface density maximum at r = 6, the 
creation of the surface density deficit for r £ [7, 7.5] and the 
new surface density maximum at r = 8.5 could be induced, 
in a similar manner as it would be by a planet, through 
the outward transport of angular momentum by the den- 
sity waves launched by the vortices. The material removed 
from the region of the deficit, which has gained angular mo- 
mentum, ends up contributing to the new surface density 
maximum at r — 8.5. 



5.7 Anticyclonic vortices 

For the disc with Q Q — 3, a vortex-pair forms at t ~ 230Po 
and lasts until the end of the simulation (t c± 300Po) and 
is reminiscent of co-orbital planets. The pair survives on 
a time-scale beyond which merging occurs in the simula- 



Figure 14. The gap profile for the Q = 3 disc at different 
times (solid, dotted and dashed) as shown through relative surface 
density perturbation. A snapshot of the profile for the disc with 
Qo = 4 (dash-dot) is also given for comparison. 



Two blobs can be identified along the gap edge, the upper 
vortex being more over-dense than the lower one. They lie 
at a radius r v = 6.4, corresponding to local surface den- 
sity maximum in azimuthally averaged one-dimensional pro- 
files, which is expected to be neutral for vortex migration 
(jPaardekooper et alJlioioh . 

is different to the pre- 



The upper vortex in Fig. 15 (a 



merging vortices or those with weak self-gravity. It has two 
spiral disturbances extending from the vortex to (r, 4>—(f P ) = 
(9, — 0.27r), whereas the pre-merging vortices have one trail- 
ing spiral. Its vortensity field is shown in Fig. |15(b)| The 
vortex core has 77 < 0, whereas the final large vortex in the 
disc with Qo = 4 has 77 > in its core, though it is still a 
local minimum. 

From the vortensity field the region where 77 < has a 
mass of M v ~ 5.46 x 10" 5 M* ~ 0.18M P and average radius 
f ~ 0.92H(r v ). This is about 18 Earth masses if M* = 1M S . 
Denoting the mean square relative velocity (to the vortex 
centre) as At; 2 , we found GM v /rAv 2 - 3.9. This region is 
gravitationally bound. Although planets of such mass are 
not expected to open gaps, the surface density deficit in 
r = [7, 7.5] in Fig. 15(a) indicates vortices may do so. This 



may be because the process is assisted by the fact that a 
vortex of size H produces a perturbation of a magnitude 
similar what would be pr oduced by Saturn when h — 0.05, 
even without self-gravity ([Paardekooper et al.ll2010h . We es- 
timate by inspection that the region with negative vorticity 
has semi-major and semi-minor axis of a ~ 1.57H(r v ), b c± 
0.55H(r v ), respectively, corresponding to aspect-ratio of 2.9. 

Finally, Fig. 15(c) show AS/S for a pair of Kida-like 
vortices placed in a disc. We note the similarity between 
the Kida-like vortex and the upper vortex in Fig. 15(a) 



tions without self- gravity. A snapshot is shown in Fig. 15 (a 



particularly the double wake structure, the tilted core and 
the surface density deficit just outside the vortices. This is a 
surprising coincidence, given that the simulation setups that 
produced them were completely different (for more details 
see below). 



6 SIMULATIONS OF CO-ORBITAL KIDA 
VORTICES 

In the simulations above, we observed one of the effects of 
self-gravity is to delay vortex merging. This effect is signif- 
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Figure 12. Vortex evolution for the disc with Q = 3. Line plot: instantaneous a viscosities (solid) and average over-density (dotted). 
Contour plots: relative surface density perturbations before and after vortex merging. 
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Figure 15. The final vortex pair at the end of the disc-planet 
simulation for Q = 3 : (a) Relative surface density contour plot 

(b) Vortensity contours for the upper vortex. Shown in (c) is the 
relative surface density for a pair of vortices formed by imposing 
Kida vortex solutions as a perturbation in a standard power-law 
disc. 



icant when vortices develop into compact structures. Fur- 
thermore, the similarity between the final vortex in the disc 
with Q = 3 and a Kida-like vortex (|Kidalll98ll ) motivates 
us to consider the effect of self-gravity on the interaction 
between two Kida vortices. 

The simulations presented here are supplementary and 
focus on the interaction between vortices without the pertur- 
bation of the planet. Thus we isolate effects due to vortex- 
vortex interactions and the influence of self-gravity on those 
without interference from the planet. In this way we obtain 
a better understanding of some of the results above. 

6.1 Numerical setup 

We use a standard power law disc with initial surface den- 
sity profile £ = Eor _1//2 , where Eo is a scaling constant 
for r £ [0.25,2.5]. The value of Eo is chosen to provide a 
specified Qxep at a specified radius as before. The initial az- 
imuthal velocity is chosen so that the disc is in hydrostatic 
equilibrium. The initial radial velocity is zero. The disc is 
locally isothermal with h — 0.05. The viscosity is v — 10 -9 
which is essentially the inviscid limit. No planets are intro- 
duced, but in this section we use (r p , ip p ) to denote a vortex's 
centroid. 

To set up v ortice s in the disc, we follow 
iLesur &; Papaloizoul (2009) and introduce velocity per- 
turbations (5u r . 5 uv) that correspond to the elliptical 
vortices of iKidal (|l98lh in an incompressible shear flow. 



Table 1. Parameters for two- vortex interaction simulations. The 
first column gives the nomenclature for the simulations. The ini- 
tial radial separation is Xq and Qi is the Keplerian Toomre 
parameter Qxep at unit radius. The fourth column indicates 
whether self-gravity was included. 



Case 


X /H 


Qi 


self-gravity 


Fnsg 





8.0 


NO 


Fsg 





8.0 


YES 


Ml 





15.9 


YES 


M2 





5.3 


YES 


SI 


0.1 


8.0 


YES 


S2 


0.2 


8.0 


YES 


S3 


0.3 


8.0 


YES 


S4 


0.6 


8.0 


YES 



Details of the implementation is described in Appendix iBl 
Two vortex perturbations are imposed with initial angular 
separation = ty/2. Table [1] summarises our numerical 
experiments. We consider switching self- gravity on or off, a 
range of disc masses (here specified through Qi = QKep(l)) 
and initial radial separations Xo of the two vortices (so that 
for one of the vortices, r p — >• r p + Xo). 

We use the FARGO code with resolution N r x = 
800 x 2400 giving ~ 17 grid points per scale height. The 
vortices are typicall y of that size. Da mping boundary con- 
ditions are applied (fde Val-Borroll2006h . 



6.1.1 Vortex pair interactions 

As an example of vortices formed by the above procedure, 
Fig. [16] shows the vortensity field for the case Fnsg. The 
vortex pair circulates at r ~ 0.96 and are localised in ra- 
dius and azimuth. The vortex centroids have vorticity —1.57 
(dimensionless units) in the non-rotating frame, close to 
the local Keplerian shear ( — 1.5r p 3 ^ 2 ~ —1.59). The vor- 
tices are not symmetric with respect to reflection in a co- 
orbital radius. The upper vortex has a long tail reaching 
the outer part of the lower vortex rather than its centroid. 
The wakes associated with each vortex are similar to those 
i nduced by a planet and responsible for vortex migration 
(jPaardekooper et alJlioioh . 

Focusing on one vortex, the region with negative ab- 
solute vorticity has half width in radius and azimuth of 
~ 0.6H and 1.3H respectively. This region has mass m v ~ 
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Figure 16. Vortices formed by imposing Kida vortex solutions 
as a perturbation to a global disc. The vortensity field (scaled) is 
shown. 



1.2 x 10~ 5 M*, or about 4 Earth masses if M* = 1M S . The 
estimate is not far from the assumption that the final vortex 
has size H with average density that of the unperturbed disc 
at the location it was set up, which gives m v ~ 1.6 x 10 -5 . 
If self- gravity is enabled, we can expect gravitational in- 
teraction between vortices to behave like that between two 
co-orbital planets of at least a few Earth masses. 



6.2 Simulation results 

We first compare cases Fsg and Fnsg. These are for the disc 
with Qi = 8 with and without self- gravity respectively. The 
distance between vortex centroids as a function of time is 
shown in Fig. flTT a). It takes almost 5 times as long for the 
self- gravitating vortices to merge. Vortensity evolutionary 
plots show that non self- gravitating vortices begin to merge 
at t = 28Po, but self- gravitating vortices at t = 126Po. For 
Fnsg, vortices approach each other within ~ 23Po of forma- 
tion and merge. Assuming they move relative to each other 
because of Keplerian shear, the time taken to merge im- 
plies that their radial separation at vortex formation must 
be c± 0.14H. This length-scale is resolved by about 2.5 grid 
cells indicating that a non-zero initial separation is gener- 
ated by grid effects, despite the vortex perturbations being 
imposed at the same radius. In a pair of lower resolution runs 
(jV r = N^/3 = 400), vortex merging without self-gravity oc- 
curs at t = 23Po whereas the with self-gravity the vortex pair 
survive until the end of the simulation. These are consistent 
if the initial radial separation is increased due to lowered 
resolution, giving stronger differential rotation , hence ear- 
lier merging in the non- self- gravity run and increased impact 
parameter in terms of horseshoe orbits in the self- gravitating 
run (see below). 

In the self- gravitating case Fsg, the inter-vortex dis- 
tance oscillates with period 50Po and there are two close- 
encounters before merging. This is analogous to the survival 
mode of the vortex pair towards the end of the disc-planet 
simulation for Q Q — 3. Given a maximal separation of c± 1.4 



and that the vortices are circulating near unit radius, the 
maximal angular separation is ~ tt/2. The minimal separa- 
tion during the first two close encounters is c± 0.6, or about 
12H which implies that the vortices are on tadpole orbits. 
As the minimum separation is much larger than the typical 
vortex size, merging does not occur during the first two close 
encounters. 

Next, in Fig. [TTT b) we vary the disc mass via Qi. The 
reference case has Qi = 8. Doubling Qi to Qi — 15.9 (case 
Ml) weakens self- gravity and vortices merge within few tens 
of Po. For Qi = 5.3 (case M2) the oscillation period is ~ 
40Po and maximal separation is 1.55, larger than for Qi = 8, 
as is the first minimum separation: the increased self- gravity 
has enhanced the mutual repulsion of vortices. The M2 two- 
vortex configuration lasts until the end of the simulation, but 
there is a secular decrease in the minimum separation (of 0.6 
at t = 20Po and 0.45 at t = 180Po) due to vortex migration. 
We expect merging to occur eventually. Consistent with the 
behaviour seen for gap edge vortices, increasing self- gravity 
delays merging. 

Our final experiment varies the initial radial separa- 
tion of the vortex perturbations. Results are shown in Fig. 
flTT c). Increasing Xo to 0.1H from the reference case, the 
first minimum separation decreases to 0.45 from 0.6: vor- 
tices approach each other more closely, but still repel and 
undergo co-orbital dynamics. The increased oscillation am- 
plitudes imply a larger tadpole orbit. For Xo = 0.2H and 
0.3H, vortex separation decreases more rapidly and become 
small enough for vortex merging. However, for Xo — 0.6 H, 
vortices simply circulate past by each other and no merg- 
ing occurs, despite reaching a similar minimal separation as 
when Xo = 0.2 H, 0.3 H. This implies merging requires that 
vortices should collide head on. 



6.3 Vortices as co-orbital planets 

The simulations presented above indicate that self- 
gravitating vortex pairs behave like co-orbital planets. This 
means there exists an initial radial separation (or impact pa- 
rameter) below which vortices execute horseshoe orbits rel- 
ative to each other. Analytical and numerical work indicates 
that vortices mer ge if their centroids are within d ^ 3s o f 
each other (e.g. IZabuskv et al.|[l979l : iMelander et alJll988h , 
where s is the vortex size. We cannot assume the same crit- 
ical d/s apply to our simulations because the situations are 
very different, but it is reasonable to expect merging if vor- 
tices can reach within some critical distance of one another. 

The above results above can be anticipated from ex- 
isting treatments of co-orbital dynami cs, which we dis- 



cuss h ere for completeness. The first is iMurrav &; Dermottl 
(2000) 's model of the co-orbital satellites of Saturn, the 
Janus- Epimetheus s ystem . The governing equation from 
IMurrav fe Dermottl (2000) gives a relationship between two 
configurations, the 'final' configuration (subscript /) and the 
'initial' configuration (subscript i) in the form 



H(0) = [sin(6'/2)]" 1 



4 

~~ 3 
2 cos 



q [H(90 
2 



H( 



(38) 



where Ar is the radial separation of the satellites, ro the 
average orbital radius (assumed fixed), their angular sepa- 
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Figure 17. Inter- vortex distances for vortex pair simulations: (a) With and without self-gravity for Qi =8 ( simulations Fsg and Fnsg) 
(b) For discs of different masses characterised by Qi, the Keplerian Toomre parameter at r = 1, self- gravity being included (simulations 
Fsg, Ml, and M2). (c) Vortices input with varying initial radial separation, Xo, with self-gravity included for disc models with Q Q = 8( 
simulations Fsg, S1,S2,S3, and S4). 



ration and q = M/M* , with M being the sum of the satellite 
masses. It is assumed g<l. 

Let us apply equation (|38|) . We assume zero final radial 
separation, Art = when the vortices are at their min- 
imum separation. Inserting Oi = 7r/2, An = 7.1 X 10 -3 
(corresponding to the initial conditions deduced for Fnsg), 
ro = 1 and q — 2.3 x 10 -5 (corresponding to the measured 
mass of the negative absolute vorticity region in Fsg) , equa- 
tion (f38|) gives 0f ~ 0.41 or a minimal angular separation 
of rsj SH so merging is not expected if vortices have size H. 
This estimate is lower than the observed value of 0.6, but we 
have used a lower limit on q. Inserting q — 1.2 x 10 -4 gives 
Of — 0.6 assuming all other parameters remain the same. 
This implies the gravitational mass of the vortex should be 
20 Earth masses. 

Equation ([38]) is illustrated in Fig. H5T a). For fixed Xo, 
as we increase q, (e.g. as a by product of increasing the disc 
mass) the minimal distance Y increases, eventually becom- 
ing too large for merging to occur. For fixed q, Y decreases 
as Xo increases. This is similar to a test particle on a horse- 
shoe orbit in the restricted three body problem. The larger 
its impact parameter (while still in libration), the closer it 
approaches the secondary mass. 

We may also view the horseshoe turns in a shearing 
sheet approximation. In Appendix \C\ we derive the inequal- 
ity 



x < 



1 / Sq q 2 
h \3yh 3y 4 h 4 



1/2 



(39) 



where xo = Xo/H, y = Y/H. Given Xo, equation ( [39]) can 
be inverted to give the range of possible Y. This inequality 
is displayed in Fig. H5T b). There exists a critical Xo = X s 
beyond which vortices do not execute horse-shoe turns and 
instead circulate past each other. There may also exist a 
value Xo, c such that Xo, c < X s and for Xo, c < Xo < X s , 
vortices execute U-turns but the values of Y attained are 
sufficiently small to allow merging. This happens only when 
the vortices have small enough mass. When Xo < Xo, c , Y 
can be larger than critical, so that merging does not occur. 

These simple models give qualitatively the same results 
as vortex pair and disc-planet simulations, and serve to ex- 
plain the effect of self-gravity delaying merging by interpret- 
ing vortices as co-orbital planets that execute mutual horse- 
shoe turns, thereby imposing a minimum inter- vortex sepa- 
ration (mainly in the azimuthal direction). If this minimum 
separation is still larger than a critical separation known to 
exist for vortex merging, then merging cannot occur. Hence, 
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Figure 18. Merging conditions based on Murrav & Der mottl 's 
model of co-orbital satellites (a) and shearing sheet dynamics 
(b). In (a), Y is the minimum inter- vortex separation and the 
horizontal line is a hypothetical critical separation below which 
merging occurs. Thus merging is less likely for larger q. Fig. llSf b) 
illustrates equation (|39|) : for a given Xo, the allowed values of min- 
imum separation Y lies between the intersection of the horizontal 
line Xo/H = constant and the solid curve. The vertical dashed 
line is a hypothetical critical separation, below which merging 
occurs. 



multiple- vortex configurations can be sustained longer as the 
strength of self- gravity is increased. 



7 IMPLICATIONS ON VORTEX-INDUCED 
MIGRATION 

Vortices at gap edges can lead to brief phases of rapid in- 
wards migration. The case for Saturn and Jupiter mass 
planets in non- self- gra vitating discs was investigated by 
iLin fe Papaloizoul (2010). In their simulations, gap edges be- 
come unstable to vortex modes leading to vortex mergers on 
dynamical time-scales, in turn resulting in a large-scale vor- 
tex circulating at either gap edge. Upon approaching the 
planet, the inner vortex can execute a horseshoe turn and 
move outwards across the gap, gaining angular momentum. 
Thus, the planet loses angular momentum and is scattered 
inwards. 

For comp arison purposes, we repeated 
ILin Sz Papaloizoul 's simulations with self-gravity. These 
correspond to inviscid discs with initially uniform surface 
density. Fig. [19] shows the orbital migration of a Saturn- 
mass planet in discs with total mass Md = 0.035M* and 
0.025M*. In the M d = 0.035M* disc with self-gravity 
neglected, vortex-induced migration occurs at t = 60Po 
whereas in the self- gravitating case it is delayed to t = 85Po . 
The delay is consistent with both the stabilisation of low 
m modes, and the slower vortex merging induced by self- 
gravity. Furthermore, there could be gravitational influence 
of the outer gap edge vortices on the inner gap edge. 
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Indeed, although somewhat inhibited , we found the 
formation of an azimuthally extended coherent vortex at the 
inner gap edge when self- gravity was included, just as in the 
non self- gravitating case.. However, it takes a longer time for 
the inner vortex to build up, disrupt the co-orbital region 
and flow across the gap, therefore induced rapid migration 
is delayed. For M d = 0.025M* it is delayed by ~ 50P , 
or twice of the higher mass disc, indicating stabilisation is 
more effective for the lower disc mass. However, the extent 
of rapid migration is unaffected by self-gravity. Notice also 
the increased oscillations when self- gravity is included. This 
is because of the sustained multi-vortex configuration at 
the outer gap edge (causing large oscillations in disc-planet 
torques), whereas without self-gravity these vortices would 
have merged. 

After the first scattering event, migration stalls while 
the planet opens a gap at its orbital radius and vortex for- 
mation recurs. In the non self- gravitating case, a large in- 
ner vortex develops (taking about 15Po) and each time it 
passes by the planet, some vortex material splits off and 
flows across the planet's orbital radius, leading to inward 
planet migration. With self gravity included, the formation 
of a single vortex takes longer (about 35Po in the higher 
mass disc) and it is narrower in radial extent than the non 
self- gravitating case. We found that in this case, when the 
thinner vortex passes by the planet, little vortex material 
splits off from the main vortex and flows across the gap. This 
may be due to the self- gravity of the vortex. Hence, there is 
a longer stalling period when self- gravity is included. Thus 
the net effect of self-gravity is to slow the migration in this 
example. 

For the setup of iLin fc Papaloizou! (|2010h . we do not 
observe a second fast migration episode within the simula- 
tion time-scale, but it may eventually occur. The total prac- 
tical simulation time of a few hundred orbits is still very 
short compared to disc lifetimes. However, for the Q = 4 
disc model used in previous parts of this work, two episodes 
of rapid migration occur. Self- gravity does not change the 
physical nature of vortex-induced migration, but we com- 
ment that for discs of even larger mass where vortex forma- 
tion no longer occurs, being replaced by global spiral insta- 
bilities, the migratio n may be even faster th an in the non 
self- gravitating case (|Lin fc Papaloizoull201ll ). 

A detailed numerical study of the effect of self-gravity 
on vortex-induced migration and its consequences for disc 
structure will be deferred to a future paper, but the simple 
experiments described here indicates that episodes of fast 
migration punctuated by periods of slower migration are still 
expected as a consequence of vortex-induced migration. 



8 CONCLUSIONS 

We have studied the effect of self-gravity on vortex forma- 
tion and evolution in protoplanetary discs. We specifically 
considered vortex production through instability associated 
with edges of surface density dips or gaps opened by a giant 
planet. It was shown analytically that vortex forming modes 

2 This is also seen in disc-planet simulations in previous sections, 
the inner gap edge is relatively more stable than outer gap edge 
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Figure 19. Vortex-induced migration with (solid) and without 
(dotted) self-gravity. 



are stabilised by self-gravity through its effect on the linear 
mode when the background remains fixed. This aspect has 
been confirmed by linear calculations, which also showed 
that self-gravity has a de-stabilising effect through its effect 
on the background state that is set up by the perturbing 
action of an embedded planet. Linear calculations showed 
that the vortex forming modes with the highest growth rate 
shift to higher m with increasing disc mass or equivalently 
decreasing Toomre Q value. 

We performed hydrodynamic simulations that compare 
simulations with and without self-gravity for a range of 
disc masses. It was found that more vortices form as the 
disc mass increased and minimum Q value decreased in ac- 
cordance with linear calculations. However, for sufficiently 
strong self- gravity, the vortex modes are suppressed. In this 
case global spiral modes develop instead. Self- gravity was 
found to delay vortex merging. In addition, multiple- vortex 
configurations can be sustained for longer with increasing 
disc mass, allowing vortices to evolve individually. 

The nature of post-merging vortices is also affected by 
self-gravity. With weak self-gravity ( i.e.. Md < 0.024M*), a 
single vortex, extended in azimuth forms and circulates at 
the outer gap edge. However, for Md = 0.031M* the final 
configuration is a vortex pair at the outer gap edge, each 
localised in azimuth. In this case a vortex, despite containing 
on the order of 2OM0, can be gravitationally bound. The 
effect of their gravitational influence on the rest of the disc 
is to redistribute mass radially whereas for lower disc masses, 
redistribution is restricted to the same azimuth. 

We note that the internal flow in these self- gravitating, 
localised vortices adjusts so that they are not destroyed by 
the background s hear, taking on a structure similar to that 
of a Kida vortex (Kida 1981). We found such vortices form 
in discs much less massive than that required for direct disc 
fragmentation. Classical gravitational instability would pro- 
duce clumps of much larger mass than vortices in our models 
and would have difficulty surviving against shear beyond a 
few orbits. However, higher resolution calculations in three 
dimensions that relax the assumption of a locally isother- 
mal equation of state are needed to assess the viability of 
the self- gravitational collapse of a vortex. 

We performed supplementary simulations of Kida-like 
vortices to understand the effect of self-gravity on inter- 
vortex interactions in a simpler setting where it was not 
contaminated by the presence of the embedded planet. The 
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effectiveness of merging as a function of disc mass may be 
understood from existing descriptions of co-orbital horse- 
shoe dynamics. In essence it appears that pressure forces 
do not play a role, so that in a self- gravitating disc, vortex 
pairs behave like co-orbital planets and as a consequence 
there exists a minimum inter- vortex distance. Merging is 
then avoided if this minimal separation is still larger than a 
critical separation below which vortex merging occurs. 

Finally, we briefly considered the consequence of self- 
gravity on vortex-induced migration. With self-gravity, the 
resistance to forming a single large vortex, which in non- 
self- gravitating simulations is responsible for scattering the 
planet inwards, results in vortex induced migration being de- 
layed. The vortices are less effective in scattering the planet 
because they are smaller and do not disrupt the co-orbital 
region as significantly as their non self- gravitating counter- 
parts. Thus, it can be said that in the regime of disc masses 
where vortices form and are significantly affected by self- 
gravity, vortex- induced migration is slowed down. 

8.1 Outlook and caveats 

A possible issue in this work is that the introduction of a 
planet over shor t timescales of a few orbits may favour insta- 
bility. However, de Val-Borr o et al. (2007) have considered 
non- self- gravitating discs with a pre-defined gap profile and 
found similar growth rates for unstable modes as models 
without an initial gap. In addition, our linear calculations 
based on gap profiles attained in a quasi-steady state are in 
good agreement with simulations. This indicates the insta- 
bility operates on such attained profiles and is not signifi- 
cantly affected by how the gap is formed. 

For the most part planetary migration has been ne- 
glected in this work. The role of the embedded planet was 
to create a structured background state from which vortices 
develop. Hence, we expect our findings on the effect of self- 
gravity to be applicable to other types of structured features 
in a protoplanetary disc that could support vortex forming 
instabilities. Such a possibility is the boundary region be- 
tween a dead zone and active region of the disc (|Lvra et al.l 
l2009h . 

Although migration was briefly considered, we note that 
there are many technical issues associated with a numerical 
study. These include the treatment of the Roche lobe and 
the time and the nature of the release of the planet. The 
issue of the effects of the choice of softening length and disc 
viscosity should also be investigated. A detailed parameter 
survey is beyond the scope of this paper, but we hope to 
investigate migration in more detail in the future. 

Vortices can trap dust particles due to their association 
with pressure maxima. We have seen that increasing self- 
gravity leads to vortices of stronger density contrast (Fig. 
[T3|) . We then anticipate self- gravitating vortices to be more 
effective at collecting solid particles. It would be interesting 
to consider a two-component, gas and dust disc model to 
investigate the specific role of self-gravity on dust trapping. 

Another issue is that we adopted the two-dimensional 
approximation. Clearly our studies should be extended to 
three dimensions. It is known that Kida vortices with as- 
pect ratio < 4 are strongly unstable to elliptic instability 
in three dimensional, unstr uctured non- self- gravitating discs 
(|Lesur &; Papaloizoul [2009). However, results may change 



when vortices are produced by an instability arising from a 
a background structure and self- gravity is important. These 
shoul d be investigated in t he future. We point out that re- 
cently Meheu t et al . (2010) demonstrated vortex formation 
in three-dimensional simulations without self-gravity via the 
vortex forming instability. Because the importance of self- 
gravity is sensitive to the disc thermodynamics and equation 
of state, it would be desirable to incorporate a realistic en- 
ergy equation in future studies. 
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Table Al. Relationship between the initial Keplerian Toomre 
stability parameter at the outer boundary, Q Q , at the planet's 
initial orbital radius, Q p , and the disc-to-star mass ratio, for disc 
models used in disc-planet interactions. 
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APPENDIX A: PARAMETERISATION OF DISC 
MODELS FOR THE VORTEX INSTABILITY 

Disc models used to study the vortex instability are labelled 
by Qo, which corresponds to values of Q p and disc-to-star 
mass ratio given in Table [ATI 



APPENDIX B: ARTIFICIAL VORTICES IN AN 
ACCRETION DISC 

The artificial Kida-like vortices used in 33 are setup as 
follows. Consider a small patch of the disc, whose centre 
(r p , Lp p ) rotates at angular speed Q p about the primary and 
set up local Cartesian co-ordinates x — r p ((p — ip p ), y = 
r p — r. In the (x,y) frame, there exists a Kida vortex solu- 
tion for incompressible flow, whose velocity field 
(u Xl u y ) is 

3Q P Cy 3Vt p x 
W * = 2(C^I)' ^ = -20C^I) (B1) 
inside the vortex core. The ratio of the vortex semi-major 
to semi-minor axis being £ = a/b is a free parameter. This 
velocity field is such that the vorticity uj is constant in the 
rotating frame. The elliptical boundary of the vortex is de- 
fined such that uo = -3Q P (1 + C 2 )/(2C(C - 1)) = ^ -30 p /2 
inside the boundary of the vortex and uj = —3Q p /2 out- 
side. The quantity uj v = —3Q P (1 + C)/(2C(C — 1)) is then 
the vorticity of the vortex core relative to the background. 
Being negative, this corresponds to an ant icy clonic vortex. 
In order to introduce perturbations corresponding to Kida 
vortices, we impose perturbations Su r = — v y , du^ = v x in- 
side a specified elliptical boundary with an exponential de- 
cay outside. The boundary is fixed by specifying £ = 8 and 
its semi-major axis a = H(r p ) where the reference radius is 

Try — 1. 



APPENDIX C: MUTUAL HORSESHOE TURNS 
IN THE SHEARING SHEET 

We describe the gravitational interactions between two vor- 
tices in the shearing sheet. It is assumed they behave like 
point masses and that pressure forces may be neglected. We 
indicate below why this is a reasonable assumption. 

Consider a local Cartesian co-ordinate system (x,y) 
that co-rotates with a small patch of fluid with angular 
velocity Q, about the primary, at a distance r p . We have 
x = r — r p , y = r p (cp — Qt). Let (xj,yj) denote the co- 
ordinates of the centroid of the j th vortex and rrij be its 
mass (j — 1,2). Defining X = X2 — x±, Y = y2 — yi and 
M = m<2 + mi, the equations of motion give 

GMX (ci) 



X - 2QY ■ 



3Q 2 X 



R 3 



Y + 2QX 



GMY 

R 3 ' 



(C2) 



where R 2 = X 2 + Y 2 . These equations imply the constancy 
of the Jacobi constant 



J 



(C3) 



Let the initial conditions be X = Xq, y = oo,X = 0,y = 
— 3Q,Xo/2. We assume the point of closest approach occurs 
when X — Y — 0. Equating J at the initial time and at the 
time of closest approach we obtain 

GM 



-lrtxt = \x 2 



Y 



(C4) 



at the time of closest approach. Since the vortices are then 
at minimum separation, Y > 0. The y component of the 
equation of motion then implies 



i . 2 i (gmV 

Substituting X from (|C4|) the inequality becomes 



3 ic 2 <r q 



q 

8Y 4 



(C5) 



(C6) 



at minimum separation, where Xo = Xo/r p , Y = Y/r p , 
q = M/M* and we have assumed Q 2 = GM*/r p . 

Eq. IC6l is useful for the case there vortices are just able 
to undergo U turns. For fixed q, the function 



' Y 8Y 4 



(C7) 



has a maximum value at Y = (q/2) 1 ^ 3 , corresponding to the 
maximum conceivable initial separation Xo = X s , where 



v 2/3 1/3 

X s = 2 q 1 r p 



(C8) 



If Xo > X s then equation ()C6|) cannot be satisfied and there 
can be no horseshoe turns. For initial separations Xo < X s , 
(|C6[) implies that the minimal inter-vortex distance must 
exceed q 1 ^ 3 /2 (so that / > 0). Now for sufficiently large g, 
q 1 / 3 /2 will be larger than the critical separation for merging, 
so merging is avoided during the encounter. 

It is interesting to compare equation (|C8|) 
to the estimate of the h ors eshoe half-width x s 
of IPaardekooper &; Papaloizoul (|20Q9h . They found 
x s = 1.6S(q/h) 1 ^ 2 r p based on hydrodynamic simula- 
tions for low mass planets. Equating x s and X s with 
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h = 0.05 we find q = 8.9 x 10" 5 . This should give the 
minimum q for which pressure effects could b e ignored. 
Inserting this value in lMurrav fc Dermott (2000) 's model of 
co-orbital satellites gives a minimal separation of 0.58 (see 
the estimate in ^6.3|) , close to simulation results. Hence, if 
the vortex-pair interaction is purely gravitational, a single 
vortex behaves in a similar way to ~ 15 Earth masses, i.e. 
a low-mass protoplanet. 

Considering a vortex size of order H, the vortex-to-star 
mass ratio is q ~ 7ri7 2 £/M* ~ h 3 /Q. For a self- gravitating 
vortex where Q ~ 1 we have q ~ h 3 = 1.25 X 10 -4 for /i — 
0.05, slightly exceeding the threshold value above. Hence we 
expect the pressureless treatment of self- gravitating vortex- 
vortex interactions to be acceptable for the purpose of ex- 
plaining the resisted-merging of self- gravitating vortices. 



